Intro

In this script, I take the collated stomach data set and calculate aggregates (feeding ratio, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row per prey type (not prey individual). I also select only the columns I need for model fitting, join environmental covariates and cpue covariates for cod and flounder, and lastly saduria biomass densities.

Load packages & source functions

# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(patchwork)
library(viridis)
library(sdmTMB)
library(mapplots)
library(sp)
library(raster)
library(ggeffects)
library(modelr)

# Source code for lan_lot to UTM
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/lon_lat_utm.R")

# Source code for map plots
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/map_plot.R")

# Trim the map plot to better fit stomach data
xmin2 <- 303000
xmax2 <- 916000
xrange <- xmax2 - xmin2

ymin2 <- 5980000
ymax2 <- 6580000
yrange <- ymax2 - ymin2

theme_set(theme_plot())

plot_map_fc2 <- plot_map_fc +
  theme(legend.position = "bottom") +
  xlim(xmin2*1.5, xmax2*0.9) +
  ylim(ymin2*1.025, ymax2*0.98268)

# Continuous colors
options(ggplot2.continuous.colour = "viridis")

Read stomach data

small_cod_stomach <- read_csv("data/clean/small_cod_stomach.csv") %>% 
  drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>% 
  rename(density_cod_small = scod_kg_km2,
         density_cod_large = lcod_kg_km2,
         density_flounder = fle_kg_km2) %>% 
  mutate(fyear = as.factor(year),
         fquarter = as.factor(quarter),
         ices_rect = as.factor(ices_rect),
         oxy_sc = (oxy - mean(oxy)) / sd(oxy),
         temp_sc = (temp - mean(temp)) / sd(temp),
         depth_sc = (depth - mean(depth)) / sd(depth),
         density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
         density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
         density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
         haul_id = as.factor(haul_id),
         species2 = "Small cod")
#> Rows: 1153 Columns: 32
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (24): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#> 
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#> 
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#> 
#>         converted 'haul_id' from character to factor (0 new NA)
#> 
#>         new variable 'fyear' (factor) with 6 unique values and 0% NA
#> 
#>         new variable 'fquarter' (factor) with 2 unique values and 0% NA
#> 
#>         new variable 'oxy_sc' (double) with 152 unique values and 0% NA
#> 
#>         new variable 'temp_sc' (double) with 152 unique values and 0% NA
#> 
#>         new variable 'depth_sc' (double) with 64 unique values and 0% NA
#> 
#>         new variable 'density_cod_small_sc' (double) with 158 unique values and 0% NA
#> 
#>         new variable 'density_cod_large_sc' (double) with 157 unique values and 0% NA
#> 
#>         new variable 'density_flounder_sc' (double) with 159 unique values and 0% NA
#> 
#>         new variable 'species2' (character) with one unique value and 0% NA

large_cod_stomach <- read_csv("data/clean/large_cod_stomach.csv") %>% 
  drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>% 
  rename(density_cod_small = scod_kg_km2,
         density_cod_large = lcod_kg_km2,
         density_flounder = fle_kg_km2) %>% 
  mutate(fyear = as.factor(year),
         fquarter = as.factor(quarter),
         ices_rect = as.factor(ices_rect),
         oxy_sc = (oxy - mean(oxy)) / sd(oxy),
         temp_sc = (temp - mean(temp)) / sd(temp),
         depth_sc = (depth - mean(depth)) / sd(depth),
         density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
         density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
         density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
         haul_id = as.factor(haul_id),
         species2 = "Large cod")
#> Rows: 2153 Columns: 32
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (24): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#> 
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#> 
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#> 
#>         converted 'haul_id' from character to factor (0 new NA)
#> 
#>         new variable 'fyear' (factor) with 6 unique values and 0% NA
#> 
#>         new variable 'fquarter' (factor) with 2 unique values and 0% NA
#> 
#>         new variable 'oxy_sc' (double) with 164 unique values and 0% NA
#> 
#>         new variable 'temp_sc' (double) with 164 unique values and 0% NA
#> 
#>         new variable 'depth_sc' (double) with 67 unique values and 0% NA
#> 
#>         new variable 'density_cod_small_sc' (double) with 158 unique values and 0% NA
#> 
#>         new variable 'density_cod_large_sc' (double) with 171 unique values and 0% NA
#> 
#>         new variable 'density_flounder_sc' (double) with 171 unique values and 0% NA
#> 
#>         new variable 'species2' (character) with one unique value and 0% NA

flounder_stomach <- read_csv("data/clean/flounder_stomach.csv") %>% 
  drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>% 
  rename(density_cod_small = scod_kg_km2,
         density_cod_large = lcod_kg_km2,
         density_flounder = fle_kg_km2) %>% 
  mutate(fyear = as.factor(year),
         fquarter = as.factor(quarter),
         ices_rect = as.factor(ices_rect),
         oxy_sc = (oxy - mean(oxy)) / sd(oxy),
         temp_sc = (temp - mean(temp)) / sd(temp),
         depth_sc = (depth - mean(depth)) / sd(depth),
         density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
         density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
         density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
         haul_id = as.factor(haul_id),
         species2 = "Flounder")
#> Rows: 2579 Columns: 32
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (24): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#> 
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#> 
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#> 
#>         converted 'haul_id' from character to factor (0 new NA)
#> 
#>         new variable 'fyear' (factor) with 6 unique values and 0% NA
#> 
#>         new variable 'fquarter' (factor) with 2 unique values and 0% NA
#> 
#>         new variable 'oxy_sc' (double) with 153 unique values and 0% NA
#> 
#>         new variable 'temp_sc' (double) with 153 unique values and 0% NA
#> 
#>         new variable 'depth_sc' (double) with 64 unique values and 0% NA
#> 
#>         new variable 'density_cod_small_sc' (double) with 148 unique values and 0% NA
#> 
#>         new variable 'density_cod_large_sc' (double) with 160 unique values and 0% NA
#> 
#>         new variable 'density_flounder_sc' (double) with 161 unique values and 0% NA
#> 
#>         new variable 'species2' (character) with one unique value and 0% NA

# Load cache
# qwraps2::lazyload_cache_dir(path = "R/main_analysis/main_diet_analysis_cache/html")

Read the prediction grid (scales variables wrt data!)

# I need predicted cod and flounder densities for this

pred_grid_1 <- read_csv("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/data/clean/pred_grid_(1_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid_2 <- read_csv("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/data/clean/pred_grid_(2_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

pred_grid <- bind_rows(pred_grid_1, pred_grid_2)

# # Flounder
# pred_grid_fle <- pred_grid %>% 
#   filter(year %in% unique(flounder_stomach$year) & ices_rect %in% unique(flounder_stomach$ices_stat_rec)) %>% 
#   # scale pred_grid wrt data
#   mutate(fyear = as.factor(year),
#          temp_sc = (temp - mean(flounder_stomach$temp)) / sd(flounder_stomach$temp), 
#          oxy_sc = (oxy - mean(flounder_stomach$oxy)) / sd(flounder_stomach$oxy), 
#          depth_sc = (depth - mean(flounder_stomach$depth)) / sd(flounder_stomach$depth), 
#          density_cod_large_sc = (density_cod_large - mean(flounder_stomach$density_cod_large)) / sd(flounder_stomach$density_cod_large), 
#          density_cod_small_sc = (density_cod_small - mean(flounder_stomach$density_cod_small)) / sd(flounder_stomach$density_cod_small), 
#          density_flounder_sc = (density_fle - mean(flounder_stomach$density_flounder)) / sd(flounder_stomach$density_flounder))
#     
# plot_map_fc +
#   geom_point(data = pred_grid_fle, aes(X, Y)) + 
#   facet_wrap(~year)
# 
# # Cod
# pred_grid_large_cod <- pred_grid %>% 
#   filter(year %in% unique(large_cod_stomach$year) & ices_rect %in% unique(large_cod_stomach$ices_stat_rec)) %>% 
#   mutate(fyear = as.factor(year),
#          temp_sc = (temp - mean(large_cod_stomach$temp)) / sd(large_cod_stomach$temp), 
#          oxy_sc = (oxy - mean(large_cod_stomach$oxy)) / sd(large_cod_stomach$oxy), 
#          depth_sc = (depth - mean(large_cod_stomach$depth)) / sd(large_cod_stomach$depth), 
#          density_cod_large_sc = (density_cod_large - mean(large_cod_stomach$density_cod_large)) / sd(large_cod_stomach$density_cod_large), 
#          density_cod_small_sc = (density_cod_small - mean(large_cod_stomach$density_cod_small)) / sd(large_cod_stomach$density_cod_small), 
#          density_flounder_sc = (density_fle - mean(large_cod_stomach$density_flounder)) / sd(large_cod_stomach$density_flounder))
# 
# plot_map_fc +
#   geom_point(data = pred_grid_large_cod, aes(X, Y)) + 
#   facet_wrap(~year)
# 
# 
# pred_grid_small_cod <- pred_grid %>% 
#   filter(year %in% unique(small_cod_stomach$year) & ices_rect %in% unique(small_cod_stomach$ices_stat_rec)) %>% 
#   mutate(fyear = as.factor(year),
#          temp_sc = (temp - mean(small_cod_stomach$temp)) / sd(small_cod_stomach$temp), 
#          oxy_sc = (oxy - mean(small_cod_stomach$oxy)) / sd(small_cod_stomach$oxy), 
#          depth_sc = (depth - mean(small_cod_stomach$depth)) / sd(small_cod_stomach$depth), 
#          density_cod_large_sc = (density_cod_large - mean(small_cod_stomach$density_cod_large)) / sd(small_cod_stomach$density_cod_large), 
#          density_cod_small_sc = (density_cod_small - mean(small_cod_stomach$density_cod_small)) / sd(small_cod_stomach$density_cod_small), 
#          density_flounder_sc = (density_fle - mean(small_cod_stomach$density_flounder)) / sd(small_cod_stomach$density_flounder))
# 
# plot_map_fc +
#   geom_point(data = pred_grid_small_cod, aes(X, Y)) + 
#   facet_wrap(~year)

Plot data in space

full_dat <- bind_rows(small_cod_stomach, large_cod_stomach, flounder_stomach) %>% 
  pivot_longer(c(tot_feeding_ratio, saduria_feeding_ratio, benthic_feeding_ratio),
               names_to = "prey_group",
               values_to = "feeding_ratio")
#> pivot_longer: reorganized (tot_feeding_ratio, benthic_feeding_ratio, saduria_feeding_ratio) into (prey_group, feeding_ratio) [was 5885x41, now 17655x40]

## Total feeding ratio
# Q1
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 1 & prey_group == "tot_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt") +
  labs(title = "Total feeding ratio", subtitle = "Quarter 1") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining


# Q4
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 4 & prey_group == "tot_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt") +
  labs(title = "Total feeding ratio", subtitle = "Quarter 4") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining


## Benthic feeding ratio
# Q1
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 1 & prey_group == "benthic_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt") +
  labs(title = "Benthic feeding ratio", subtitle = "Quarter 1") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining


# Q4
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 4 & prey_group == "benthic_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt") +
  labs(title = "Benthic feeding ratio", subtitle = "Quarter 4") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining


## Saduria feeding ratio
# Q1
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 1 & prey_group == "saduria_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt") +
  labs(title = "Saduria feeding ratio", subtitle = "Quarter 1") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining


# Q4
plot_map_fc2 +
  geom_point(data = filter(full_dat, quarter == 4 & prey_group == "saduria_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
  facet_grid(species2 ~ year) +
  scale_color_viridis(trans = "sqrt") +
  labs(title = "Saduria feeding ratio", subtitle = "Quarter 4") +
  theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining

Plot data (response vs explanatory variables)

ggplot(full_dat, aes(depth, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(full_dat, aes(oxy, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(full_dat, aes(temp, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(full_dat, aes(density_flounder_sc, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(full_dat, aes(density_cod_large_sc, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'


ggplot(full_dat, aes(density_cod_small_sc, feeding_ratio)) + 
  geom_point() + 
  stat_smooth() + 
  facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Fit spatiotemporal model to feeding ratio

sdmTMB models with densities as covariates to stomach content data. First make mesh:

# Small cod
small_cod_spde <- make_mesh(small_cod_stomach,
                            c("X", "Y"),
                            cutoff = 15,
                            type = "kmeans",
                            seed = 42)
plot(small_cod_spde)



# Large cod
large_cod_spde <- make_mesh(large_cod_stomach,
                            c("X", "Y"),
                            n_knots = 75,
                            type = "kmeans",
                            seed = 42)
plot(large_cod_spde)



# Flounder
flounder_spde <- make_mesh(flounder_stomach,
                           c("X", "Y"),
                           n_knots = 75,
                           type = "kmeans",
                           seed = 42)
plot(flounder_spde)

# s_cod_tot0 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter,
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "off",
#   mesh = small_cod_spde,
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# tidy(s_cod_tot0)
# sanity(s_cod_tot0)
# 
# # Plot residuals in space
# small_cod_stomach$tot_resid0 <- residuals(s_cod_tot0)
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid0), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# 
# s_cod_tot <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "off",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# tidy(s_cod_tot)
# sanity(s_cod_tot)
# 
# # Plot residuals in space
# small_cod_stomach$tot_resid <- residuals(s_cod_tot)
# 
# # plot_map_fc2 +
# #   geom_point(data = small_cod_stomach, aes(X*1000, Y*1000, fill = tot_resid), size = 3,
# #              shape = 21, color = "grey50", stroke = 0.2) +
# #   facet_wrap(~ year) +
# #   scale_fill_gradient2()
# #
# # plot_map_fc2 +
# #   geom_point(data = filter(small_cod_stomach, tot_resid < 5), aes(X*1000, Y*1000, fill = tot_resid), size = 3,
# #              shape = 21, color = "grey50", stroke = 0.2) +
# #   facet_wrap(~ year) +
# #   scale_fill_gradient2()
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# # Fit model now with spatial random effect
# s_cod_tot2 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "on",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# print(s_cod_tot2)
# tidy(s_cod_tot2)
# sanity(s_cod_tot2)
# 
# small_cod_stomach$tot_resid2 <- residuals(s_cod_tot2)
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid2), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# # Fit model now with spatiotemporal random effect
# s_cod_tot3 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "IID",
#   spatial = "off",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# print(s_cod_tot3)
# tidy(s_cod_tot3)
# sanity(s_cod_tot3)
# 
# small_cod_stomach$tot_resid3 <- residuals(s_cod_tot3)
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid3), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# # Fit model now with spatiotemporal AND spatial random effect
# s_cod_tot4 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "IID",
#   spatial = "on",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# print(s_cod_tot4)
# tidy(s_cod_tot4)
# sanity(s_cod_tot4)
# 
# small_cod_stomach$tot_resid4 <- residuals(s_cod_tot4)
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid4), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# 
# # Just random effects...
# small_cod_stomach$ices_rect <- as.factor(small_cod_stomach$ices_rect)
# 
# s_cod_tot5 <- sdmTMB(
#   data = small_cod_stomach,
#   formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|ices_rect) + (1|haul_id),
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "off",
#   mesh = small_cod_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )
# 
# print(s_cod_tot5)
# tidy(s_cod_tot5)
# sanity(s_cod_tot5)
# 
# small_cod_stomach$tot_resid5 <- residuals(s_cod_tot5)
# 
# plot_map_fc2 +
#   geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid5), size = 3,
#              shape = 21, color = "grey50", stroke = 0.2) +
#   scale_fill_gradient2()
# 
# 
# # Plot together
# long <- small_cod_stomach %>%
#   pivot_longer(c(tot_resid0, tot_resid, tot_resid2, tot_resid3, tot_resid4, tot_resid5), names_to = "Model", values_to = "Residuals")
# 
# plot_map_fc2 +
#   geom_jitter(data = long %>% filter(Residuals < 4), aes(X*1000, Y*1000, fill = Residuals), size = 3,
#              shape = 21, color = "grey50", stroke = 0.1) +
#   scale_fill_gradient2() +
#   facet_wrap(~Model, ncol = 3)
# 
# plot_map_fc2 +
#   geom_jitter(data = long %>% filter(Residuals < 4), aes(X*1000, Y*1000, fill = Residuals), size = 3,
#              shape = 21, color = "grey50", stroke = 0.1) +
#   scale_fill_gradient2() +
#   facet_grid(year ~ Model)
# 
# 
# # Plot spatial and spatiotemporal random effect
# pred_grid_small <- pred_grid %>%
#   filter(quarter == 4) %>%
#   mutate(X = X/1000,
#          Y = Y/1000) %>%
#   filter(year %in% unique(small_cod_stomach$year)) %>%
#   filter(X < max(small_cod_stomach$X)) %>%
#   filter(X > min(small_cod_stomach$X)) %>%
#   filter(Y < max(small_cod_stomach$Y)) %>%
#   filter(Y > min(small_cod_stomach$Y))
# 
# pred_grid_small <- pred_grid_small %>%
#   mutate(fyear = factor(year),
#          fquarter = factor(quarter),
#          temp_sc = 0,
#          depth_sc = 0,
#          density_cod_small_sc = 0,
#          density_cod_large_sc = 0,
#          density_flounder_sc = 0,
#          oxy_sc = 0)
# 
# spat_pred <- predict(s_cod_tot2, newdata = pred_grid_small)
# 
# plot_map_fc2 +
#   geom_raster(data = spat_pred, aes(X*1000, Y*1000, fill = omega_s)) +
#   scale_fill_gradient2()
# 
# spat_pred2 <- predict(s_cod_tot3, newdata = pred_grid_small)
# 
# plot_map_fc2 +
#   geom_raster(data = spat_pred2, aes(X*1000, Y*1000, fill = epsilon_st)) +
#   scale_fill_gradient2() +
#   facet_wrap(~year)
# 
# AIC(s_cod_tot0, s_cod_tot, s_cod_tot2, s_cod_tot3, s_cod_tot4, s_cod_tot5)

Small cod

# Total
s_cod_tot <- sdmTMB(
  data = small_cod_stomach,
  formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = small_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(s_cod_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     density_cod_small_sc + density_cod_large_sc + density_flounder_sc * 
#>  Formula:     oxy_sc + (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: tweedie(link = 'log')
#>  
#>                            coef.est coef.se
#> fyear2015                     -5.27    0.41
#> fyear2016                     -5.42    0.26
#> fyear2017                     -4.72    0.20
#> fyear2018                     -5.57    0.21
#> fyear2019                     -4.87    0.37
#> fyear2020                     -5.59    0.19
#> fquarter4                     -0.09    0.31
#> temp_sc                       -0.18    0.11
#> depth_sc                      -0.31    0.14
#> density_cod_small_sc          -0.11    0.23
#> density_cod_large_sc           0.17    0.20
#> density_flounder_sc            0.04    0.10
#> oxy_sc                        -0.12    0.16
#> density_flounder_sc:oxy_sc    -0.22    0.12
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.77
#> 
#> Dispersion parameter: 0.18
#> Tweedie p: 1.59
#> ML criterion at convergence: -2952.284
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_tot)
#>                          term    estimate  std.error
#> 1                   fyear2015 -5.27421186 0.41041126
#> 2                   fyear2016 -5.42038244 0.26174994
#> 3                   fyear2017 -4.71694581 0.19506241
#> 4                   fyear2018 -5.56547652 0.20792756
#> 5                   fyear2019 -4.87241401 0.37423092
#> 6                   fyear2020 -5.58949458 0.18710427
#> 7                   fquarter4 -0.09295426 0.31180598
#> 8                     temp_sc -0.18332144 0.11378328
#> 9                    depth_sc -0.31052577 0.14280439
#> 10       density_cod_small_sc -0.11053333 0.22523374
#> 11       density_cod_large_sc  0.16726527 0.19789911
#> 12        density_flounder_sc  0.04175951 0.09689319
#> 13                     oxy_sc -0.12249743 0.15996736
#> 14 density_flounder_sc:oxy_sc -0.22168059 0.11837020
sanity(s_cod_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Benthic
s_cod_ben <- sdmTMB(
  data = small_cod_stomach,
  formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = small_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(s_cod_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     density_cod_small_sc + density_cod_large_sc + density_flounder_sc * 
#>  Formula:     oxy_sc + (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: tweedie(link = 'log')
#>  
#>                            coef.est coef.se
#> fyear2015                     -5.64    0.42
#> fyear2016                     -5.72    0.27
#> fyear2017                     -4.97    0.20
#> fyear2018                     -5.67    0.21
#> fyear2019                     -5.55    0.39
#> fyear2020                     -5.90    0.19
#> fquarter4                      0.22    0.32
#> temp_sc                       -0.18    0.12
#> depth_sc                      -0.36    0.15
#> density_cod_small_sc           0.13    0.23
#> density_cod_large_sc          -0.13    0.21
#> density_flounder_sc            0.10    0.10
#> oxy_sc                         0.10    0.17
#> density_flounder_sc:oxy_sc    -0.29    0.12
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.79
#> 
#> Dispersion parameter: 0.14
#> Tweedie p: 1.56
#> ML criterion at convergence: -2943.031
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_ben)
#>                          term    estimate  std.error
#> 1                   fyear2015 -5.63963437 0.41728201
#> 2                   fyear2016 -5.71555073 0.26736376
#> 3                   fyear2017 -4.96931045 0.20133219
#> 4                   fyear2018 -5.67299371 0.21258020
#> 5                   fyear2019 -5.54980300 0.38532858
#> 6                   fyear2020 -5.89510007 0.19336637
#> 7                   fquarter4  0.21970967 0.32106970
#> 8                     temp_sc -0.17986537 0.11810696
#> 9                    depth_sc -0.35792181 0.14677782
#> 10       density_cod_small_sc  0.12910755 0.22899007
#> 11       density_cod_large_sc -0.12685493 0.20658268
#> 12        density_flounder_sc  0.09636279 0.09872993
#> 13                     oxy_sc  0.09878458 0.16518621
#> 14 density_flounder_sc:oxy_sc -0.29431840 0.12153496
sanity(s_cod_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Saduria
s_cod_sad <- sdmTMB(
  data = small_cod_stomach,
  formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = small_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(s_cod_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     density_cod_small_sc + density_cod_large_sc + density_flounder_sc * 
#>  Formula:     oxy_sc + (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: tweedie(link = 'log')
#>  
#>                            coef.est coef.se
#> fyear2015                     -9.08    1.29
#> fyear2016                     -8.99    0.84
#> fyear2017                    -10.53    0.74
#> fyear2018                    -10.22    0.79
#> fyear2019                    -10.75    1.33
#> fyear2020                     -9.04    0.59
#> fquarter4                      0.34    1.06
#> temp_sc                        0.14    0.33
#> depth_sc                       0.26    0.40
#> density_cod_small_sc          -0.15    0.54
#> density_cod_large_sc           0.08    0.55
#> density_flounder_sc           -1.34    0.45
#> oxy_sc                         0.35    0.46
#> density_flounder_sc:oxy_sc    -0.41    0.49
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      2.01
#> 
#> Dispersion parameter: 0.19
#> Tweedie p: 1.48
#> ML criterion at convergence: -263.206
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_sad)
#>                          term     estimate std.error
#> 1                   fyear2015  -9.07779271 1.2898235
#> 2                   fyear2016  -8.98853609 0.8428161
#> 3                   fyear2017 -10.52518183 0.7365606
#> 4                   fyear2018 -10.21759621 0.7939155
#> 5                   fyear2019 -10.75111413 1.3345506
#> 6                   fyear2020  -9.04136144 0.5937819
#> 7                   fquarter4   0.34186222 1.0616702
#> 8                     temp_sc   0.14426222 0.3303548
#> 9                    depth_sc   0.26102087 0.4041756
#> 10       density_cod_small_sc  -0.15380459 0.5352128
#> 11       density_cod_large_sc   0.08468962 0.5467832
#> 12        density_flounder_sc  -1.34276867 0.4522182
#> 13                     oxy_sc   0.35260152 0.4586600
#> 14 density_flounder_sc:oxy_sc  -0.40751860 0.4911995
sanity(s_cod_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

Large cod

# Total
l_cod_tot <- sdmTMB(
  data = large_cod_stomach,
  formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = large_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(l_cod_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     density_cod_small_sc + density_cod_large_sc + density_flounder_sc * 
#>  Formula:     oxy_sc + (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: tweedie(link = 'log')
#>  
#>                            coef.est coef.se
#> fyear2015                     -3.93    0.24
#> fyear2016                     -4.49    0.16
#> fyear2017                     -4.25    0.13
#> fyear2018                     -4.33    0.11
#> fyear2019                     -4.44    0.24
#> fyear2020                     -4.56    0.13
#> fquarter4                     -0.48    0.20
#> temp_sc                        0.04    0.08
#> depth_sc                       0.05    0.11
#> density_cod_small_sc           0.42    0.22
#> density_cod_large_sc          -0.22    0.22
#> density_flounder_sc           -0.06    0.06
#> oxy_sc                         0.06    0.10
#> density_flounder_sc:oxy_sc    -0.03    0.07
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.42
#> 
#> Dispersion parameter: 0.60
#> Tweedie p: 1.71
#> ML criterion at convergence: -4812.982
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_tot)
#>                          term    estimate  std.error
#> 1                   fyear2015 -3.93309844 0.24356798
#> 2                   fyear2016 -4.49285009 0.15549512
#> 3                   fyear2017 -4.25319947 0.13076704
#> 4                   fyear2018 -4.32775670 0.11248077
#> 5                   fyear2019 -4.44268256 0.24316977
#> 6                   fyear2020 -4.56147977 0.12705887
#> 7                   fquarter4 -0.47584719 0.20471977
#> 8                     temp_sc  0.03629594 0.07886520
#> 9                    depth_sc  0.05126377 0.10663245
#> 10       density_cod_small_sc  0.41917028 0.21975227
#> 11       density_cod_large_sc -0.22222104 0.21800609
#> 12        density_flounder_sc -0.05704968 0.06125595
#> 13                     oxy_sc  0.06168315 0.10263081
#> 14 density_flounder_sc:oxy_sc -0.03314841 0.06987393
sanity(l_cod_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Benthic
l_cod_ben <- sdmTMB(
  data = large_cod_stomach,
  formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = large_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(l_cod_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     density_cod_small_sc + density_cod_large_sc + density_flounder_sc * 
#>  Formula:     oxy_sc + (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: tweedie(link = 'log')
#>  
#>                            coef.est coef.se
#> fyear2015                     -6.81    0.39
#> fyear2016                     -6.65    0.25
#> fyear2017                     -6.55    0.21
#> fyear2018                     -6.66    0.18
#> fyear2019                     -7.78    0.38
#> fyear2020                     -6.79    0.20
#> fquarter4                      1.07    0.32
#> temp_sc                       -0.23    0.12
#> depth_sc                       0.15    0.17
#> density_cod_small_sc           0.05    0.33
#> density_cod_large_sc           0.29    0.33
#> density_flounder_sc           -0.14    0.10
#> oxy_sc                         0.34    0.17
#> density_flounder_sc:oxy_sc    -0.01    0.12
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.79
#> 
#> Dispersion parameter: 0.42
#> Tweedie p: 1.67
#> ML criterion at convergence: -4579.306
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_ben)
#>                          term     estimate  std.error
#> 1                   fyear2015 -6.810709941 0.39070269
#> 2                   fyear2016 -6.650622983 0.25055841
#> 3                   fyear2017 -6.548153215 0.20997551
#> 4                   fyear2018 -6.659105762 0.18475239
#> 5                   fyear2019 -7.778384365 0.38493973
#> 6                   fyear2020 -6.786058885 0.20188801
#> 7                   fquarter4  1.071428531 0.32284539
#> 8                     temp_sc -0.230415833 0.12383263
#> 9                    depth_sc  0.149296568 0.17155416
#> 10       density_cod_small_sc  0.047376615 0.33401962
#> 11       density_cod_large_sc  0.285332497 0.33497684
#> 12        density_flounder_sc -0.135768118 0.09772925
#> 13                     oxy_sc  0.343186202 0.16712306
#> 14 density_flounder_sc:oxy_sc -0.006315457 0.11754496
sanity(l_cod_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Saduria
l_cod_sad <- sdmTMB(
  data = large_cod_stomach,
  formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = large_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(l_cod_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     density_cod_small_sc + density_cod_large_sc + density_flounder_sc * 
#>  Formula:     oxy_sc + (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: tweedie(link = 'log')
#>  
#>                            coef.est coef.se
#> fyear2015                     -8.84    1.26
#> fyear2016                    -10.91    0.80
#> fyear2017                    -11.11    0.75
#> fyear2018                    -11.28    0.67
#> fyear2019                    -12.09    1.40
#> fyear2020                    -10.54    0.64
#> fquarter4                     -0.48    1.04
#> temp_sc                        0.25    0.35
#> depth_sc                       0.22    0.48
#> density_cod_small_sc           0.03    0.65
#> density_cod_large_sc          -0.45    0.71
#> density_flounder_sc           -1.08    0.41
#> oxy_sc                         0.68    0.47
#> density_flounder_sc:oxy_sc     0.00    0.47
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      2.45
#> 
#> Dispersion parameter: 0.12
#> Tweedie p: 1.46
#> ML criterion at convergence: -487.780
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_sad)
#>                          term      estimate std.error
#> 1                   fyear2015  -8.843930672 1.2551682
#> 2                   fyear2016 -10.910267097 0.8005328
#> 3                   fyear2017 -11.107997362 0.7513374
#> 4                   fyear2018 -11.278581130 0.6708063
#> 5                   fyear2019 -12.089648568 1.4039963
#> 6                   fyear2020 -10.540081776 0.6418378
#> 7                   fquarter4  -0.477676767 1.0426728
#> 8                     temp_sc   0.253550124 0.3549809
#> 9                    depth_sc   0.218380287 0.4845376
#> 10       density_cod_small_sc   0.030106804 0.6508573
#> 11       density_cod_large_sc  -0.451303087 0.7126573
#> 12        density_flounder_sc  -1.083810480 0.4104634
#> 13                     oxy_sc   0.675276134 0.4674646
#> 14 density_flounder_sc:oxy_sc   0.004382066 0.4720712
sanity(l_cod_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

Flounder

# Total
fle_tot <- sdmTMB(
  data = flounder_stomach,
  formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = flounder_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(fle_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     density_cod_small_sc * oxy_sc + density_cod_large_sc + density_flounder_sc + 
#>  Formula:     (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: tweedie(link = 'log')
#>  
#>                             coef.est coef.se
#> fyear2015                      -5.07    0.33
#> fyear2016                      -4.96    0.26
#> fyear2017                      -4.95    0.17
#> fyear2018                      -5.16    0.17
#> fyear2019                      -5.01    0.30
#> fyear2020                      -4.39    0.14
#> fquarter4                       0.47    0.25
#> temp_sc                         0.05    0.09
#> depth_sc                       -0.13    0.12
#> density_cod_small_sc            0.27    0.26
#> oxy_sc                          0.21    0.14
#> density_cod_large_sc           -0.25    0.26
#> density_flounder_sc            -0.08    0.08
#> density_cod_small_sc:oxy_sc    -0.13    0.17
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.69
#> 
#> Dispersion parameter: 0.17
#> Tweedie p: 1.49
#> ML criterion at convergence: -3635.533
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_tot)
#>                           term    estimate  std.error
#> 1                    fyear2015 -5.06715201 0.32774948
#> 2                    fyear2016 -4.96311363 0.25612380
#> 3                    fyear2017 -4.94872102 0.16634274
#> 4                    fyear2018 -5.15532320 0.17090355
#> 5                    fyear2019 -5.01188780 0.30354649
#> 6                    fyear2020 -4.39318789 0.13752021
#> 7                    fquarter4  0.47129234 0.25039002
#> 8                      temp_sc  0.05425746 0.09198798
#> 9                     depth_sc -0.13324303 0.11759711
#> 10        density_cod_small_sc  0.26816835 0.26459501
#> 11                      oxy_sc  0.21363168 0.13699214
#> 12        density_cod_large_sc -0.25338400 0.26356414
#> 13         density_flounder_sc -0.08200934 0.07529710
#> 14 density_cod_small_sc:oxy_sc -0.13353163 0.16917662
sanity(fle_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Benthic
fle_ben <- sdmTMB(
  data = flounder_stomach,
  formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = flounder_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(fle_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     density_cod_small_sc * oxy_sc + density_cod_large_sc + density_flounder_sc + 
#>  Formula:     (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: tweedie(link = 'log')
#>  
#>                             coef.est coef.se
#> fyear2015                      -5.11    0.34
#> fyear2016                      -5.00    0.27
#> fyear2017                      -4.96    0.17
#> fyear2018                      -5.13    0.18
#> fyear2019                      -5.30    0.32
#> fyear2020                      -4.45    0.14
#> fquarter4                       0.47    0.26
#> temp_sc                         0.10    0.10
#> depth_sc                       -0.20    0.12
#> density_cod_small_sc            0.32    0.27
#> oxy_sc                          0.23    0.14
#> density_cod_large_sc           -0.29    0.27
#> density_flounder_sc            -0.06    0.08
#> density_cod_small_sc:oxy_sc    -0.09    0.18
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      0.72
#> 
#> Dispersion parameter: 0.17
#> Tweedie p: 1.49
#> ML criterion at convergence: -3574.297
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_ben)
#>                           term    estimate  std.error
#> 1                    fyear2015 -5.11048157 0.33989780
#> 2                    fyear2016 -5.00496148 0.26663702
#> 3                    fyear2017 -4.96251227 0.17279852
#> 4                    fyear2018 -5.12786504 0.17730557
#> 5                    fyear2019 -5.29751455 0.31926186
#> 6                    fyear2020 -4.44609046 0.14312220
#> 7                    fquarter4  0.46914758 0.26092314
#> 8                      temp_sc  0.09771251 0.09572682
#> 9                     depth_sc -0.19742836 0.12302319
#> 10        density_cod_small_sc  0.31634925 0.27402938
#> 11                      oxy_sc  0.22536646 0.14331473
#> 12        density_cod_large_sc -0.28796091 0.27286263
#> 13         density_flounder_sc -0.05528865 0.07821425
#> 14 density_cod_small_sc:oxy_sc -0.08752974 0.17644492
sanity(fle_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100

# Saduria
fle_sad <- sdmTMB(
  data = flounder_stomach,
  formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = flounder_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

print(fle_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + 
#>  Formula:     density_cod_small_sc * oxy_sc + density_cod_large_sc + density_flounder_sc + 
#>  Formula:     (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: tweedie(link = 'log')
#>  
#>                             coef.est coef.se
#> fyear2015                      -6.33    1.21
#> fyear2016                      -5.62    0.91
#> fyear2017                      -8.61    0.61
#> fyear2018                     -10.12    0.69
#> fyear2019                      -8.55    1.18
#> fyear2020                      -8.48    0.56
#> fquarter4                      -1.26    0.91
#> temp_sc                        -0.40    0.31
#> depth_sc                        0.03    0.37
#> density_cod_small_sc           -0.98    0.63
#> oxy_sc                          0.26    0.43
#> density_cod_large_sc            0.51    0.62
#> density_flounder_sc            -1.00    0.37
#> density_cod_small_sc:oxy_sc    -0.35    0.56
#> 
#> Random intercepts:
#>         Std. Dev.
#> haul_id      2.63
#> 
#> Dispersion parameter: 0.21
#> Tweedie p: 1.49
#> ML criterion at convergence: -1086.186
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_sad)
#>                           term     estimate std.error
#> 1                    fyear2015  -6.32974501 1.2138097
#> 2                    fyear2016  -5.62234282 0.9086070
#> 3                    fyear2017  -8.60780551 0.6072356
#> 4                    fyear2018 -10.12096544 0.6889580
#> 5                    fyear2019  -8.55498870 1.1832802
#> 6                    fyear2020  -8.48135360 0.5636020
#> 7                    fquarter4  -1.25920836 0.9118550
#> 8                      temp_sc  -0.40038843 0.3143738
#> 9                     depth_sc   0.03433966 0.3686257
#> 10        density_cod_small_sc  -0.98458983 0.6321219
#> 11                      oxy_sc   0.26364504 0.4285877
#> 12        density_cod_large_sc   0.50812859 0.6226676
#> 13         density_flounder_sc  -1.00034650 0.3727599
#> 14 density_cod_small_sc:oxy_sc  -0.35448644 0.5617143
sanity(fle_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# mcmc_res_l_cod_sad_v2 <- residuals(l_cod_sad_v2,
#                                    type = "mle-mcmc",
#                                    mcmc_iter = 201,
#                                    mcmc_warmup = 200)
# 
# qqnorm(mcmc_res_l_cod_sad_v2, ylim = c(-5, 5), xlim = c(-5, 5)); qqline(mcmc_res_l_cod_sad_v2)

Plot coefficients

#### Total prey
# small cod total prey
small_cod_tot_pars <- bind_rows(tidy(s_cod_tot, effects = "ran_pars", conf.int = TRUE),
                                tidy(s_cod_tot, effects = "fixed", conf.int = TRUE)) %>%
  mutate(species = "Small cod",
         response = "Total prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# large cod total prey
large_cod_tot_pars <- bind_rows(tidy(l_cod_tot, effects = "ran_pars", conf.int = TRUE),
                                tidy(l_cod_tot, effects = "fixed", conf.int = TRUE)) %>%
  mutate(species = "Large cod",
         response = "Total prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# flounder total prey
flounder_tot_pars <- bind_rows(tidy(fle_tot, effects = "ran_pars", conf.int = TRUE),
                               tidy(fle_tot, effects = "fixed", conf.int = TRUE)) %>%
  mutate(species = "Flounder",
         response = "Total prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

#### Benthic prey
# small cod benthic prey
small_cod_ben_pars <- bind_rows(tidy(s_cod_ben, effects = "ran_pars", conf.int = TRUE),
                                tidy(s_cod_ben, effects = "fixed", conf.int = TRUE)) %>%
  mutate(species = "Small cod",
         response = "Benthic prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# large cod benthic prey
large_cod_ben_pars <- bind_rows(tidy(l_cod_ben, effects = "ran_pars", conf.int = TRUE),
                                tidy(l_cod_ben, effects = "fixed", conf.int = TRUE)) %>%
  mutate(species = "Large cod",
         response = "Benthic prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# flounder benthic prey
flounder_ben_pars <- bind_rows(tidy(fle_ben, effects = "ran_pars", conf.int = TRUE),
                               tidy(fle_ben, effects = "fixed", conf.int = TRUE)) %>%
  mutate(species = "Flounder",
         response = "Benthic prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA


#### Saduria
# small cod saduria
small_cod_sad_pars <- bind_rows(tidy(s_cod_sad, effects = "ran_pars", conf.int = TRUE),
                                tidy(s_cod_sad, effects = "fixed", conf.int = TRUE)) %>%
  mutate(species = "Small cod",
         response = "Saduria FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# large cod saduria
large_cod_sad_pars <- bind_rows(tidy(l_cod_sad, effects = "ran_pars", conf.int = TRUE),
                                tidy(l_cod_sad, effects = "fixed", conf.int = TRUE)) %>%
  mutate(species = "Large cod",
         response = "Saduria FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

# flounder saduria
flounder_sad_pars <- bind_rows(tidy(fle_sad, effects = "ran_pars", conf.int = TRUE),
                               tidy(fle_sad, effects = "fixed", conf.int = TRUE)) %>%
  mutate(species = "Flounder",
         response = "Saduria FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

coefs <- bind_rows(small_cod_tot_pars,
                   large_cod_tot_pars,
                   flounder_tot_pars,
                      
                   small_cod_ben_pars,
                   large_cod_ben_pars,
                   flounder_ben_pars,
                   
                   small_cod_sad_pars,
                   large_cod_sad_pars,
                   flounder_sad_pars)

coefs <- coefs %>% 
  mutate(term2 = term,
         term2 = ifelse(term == "temp_sc", "Temperature", term2),
         term2 = ifelse(term == "oxy_sc", "Oxygen", term2),
         term2 = ifelse(term == "depth_sc", "Depth", term2),
         term2 = ifelse(term == "density_cod_large_sc", "Cod density (> 29)", term2),
         term2 = ifelse(term == "density_cod_small_sc", "Cod density (< 29)", term2),
         term2 = ifelse(term == "density_cod_small_sc:oxy_sc", "Cod density (> 29) x Oxygen", term2),
         term2 = ifelse(term == "density_flounder_sc:oxy_sc", "Flounder density x Oxygen", term2),
         term2 = ifelse(term == "density_cod_small_sc:oxy_sc", "Small cod density x Oxygen", term2),
         term2 = ifelse(term == "density_flounder_sc", "Flounder density", term2))
#> mutate: new variable 'term2' (character) with 18 unique values and 0% NA

# Fixed continuous effects:
coefs %>% 
  filter(!term %in% c("phi", "sigma_G", "tweedie_p")) %>% 
  filter(!grepl('year', term)) %>% 
  filter(!grepl('quarter', term)) %>% 
  ggplot(aes(term2, estimate, color = species)) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.2) +
  geom_point(position = position_dodge(width = 0.8)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
                position = position_dodge(width = 0.8), alpha = 0.5) +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  facet_wrap(~response) +
  labs(x = "Predictor", y = "Standardized coefficient") +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90),
        aspect.ratio = 1.5) +
  NULL 
#> filter: removed 27 rows (18%), 126 rows remaining
#> filter: removed 54 rows (43%), 72 rows remaining
#> filter: removed 9 rows (12%), 63 rows remaining


ggsave("figures/fr_coefs.pdf", width = 17, height = 17, units = "cm")

coefs %>% 
  filter(term %in% c("fyear2015", "fyear2016", "fyear2017", "fyear2018", "fyear2019", "fyear2020", "quarter")) %>% 
  mutate(term2 = ifelse(term == "fyear2015", "2015", term2),
         term2 = ifelse(term == "fyear2016", "2016", term2),
         term2 = ifelse(term == "fyear2017", "2017", term2),
         term2 = ifelse(term == "fyear2018", "2018", term2),
         term2 = ifelse(term == "fyear2019", "2019", term2),
         term2 = ifelse(term == "fyear2020", "2020", term2),
         term2 = ifelse(term == "quarter", "Quarter", term2)) %>% 
  filter(!term2 == "Quarter") %>% 
  ggplot(aes(term2, estimate, color = species)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  facet_wrap(~response, scales = "free") +
  labs(x = "Predictor", y = "Standardized coefficient") +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90),
        aspect.ratio = 1.5) +
  NULL 
#> filter: removed 99 rows (65%), 54 rows remaining
#> mutate: changed 54 values (100%) of 'term2' (0 new NA)
#> filter: no rows removed


ggsave("figures/supp/fr_coefs_yr.pdf", width = 17, height = 17, units = "cm")
# fle_sad <- sdmTMB(
#   data = flounder_stomach,
#   formula = saduria_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
#   family = tweedie(link = "log"),
#   time = "year",
#   spatiotemporal = "off",
#   spatial = "off",
#   mesh = flounder_spde,
#   priors = sdmTMBpriors(
#     b = normal(c(rep(NA, 7), rep(0, 7)),
#                c(rep(NA, 7), rep(1, 7)))),
#   control = sdmTMBcontrol(newton_loops = 1)
#   )

s_cod_sad2 <- sdmTMB(
  data = small_cod_stomach,
  formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = small_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )

l_cod_sad2 <- sdmTMB(
  data = large_cod_stomach,
  formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = large_cod_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )


fle_sad2 <- sdmTMB(
  data = flounder_stomach,
  formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc,
  family = tweedie(link = "log"),
  time = "year",
  spatiotemporal = "off",
  spatial = "off",
  mesh = flounder_spde,
  priors = sdmTMBpriors(
    b = normal(c(rep(NA, 7), rep(0, 7)),
               c(rep(NA, 7), rep(1, 7)))),
  control = sdmTMBcontrol(newton_loops = 1)
  )


quantile(flounder_stomach$density_flounder, prob = c(0.1, 0.9))
#>       10%       90% 
#>  231.5209 4716.5631

nd_fle <- data.frame(density_flounder = seq(0, 5000, length.out = 100))

nd_fle <- nd_fle %>%
  mutate(density_flounder_sc = (density_flounder - mean(flounder_stomach$density_flounder)) / sd(flounder_stomach$density_flounder),
         year = 2017L,
         fyear = as.factor(2017),
         fquarter = as.factor(1),
         temp_sc = 0,
         depth_sc = 0,
         oxy_sc = 0,
         density_cod_small_sc = 0,
         density_cod_large_sc = 0
         )
#> mutate: new variable 'density_flounder_sc' (double) with 100 unique values and 0% NA
#>         new variable 'year' (integer) with one unique value and 0% NA
#>         new variable 'fyear' (factor) with one unique value and 0% NA
#>         new variable 'fquarter' (factor) with one unique value and 0% NA
#>         new variable 'temp_sc' (double) with one unique value and 0% NA
#>         new variable 'depth_sc' (double) with one unique value and 0% NA
#>         new variable 'oxy_sc' (double) with one unique value and 0% NA
#>         new variable 'density_cod_small_sc' (double) with one unique value and 0% NA
#>         new variable 'density_cod_large_sc' (double) with one unique value and 0% NA

# Predict from full model
sc_margin_sad <- predict(s_cod_sad2, newdata = nd_fle, se_fit = TRUE, re_form = NA, re_form_iid = NA) %>% 
  mutate(species2 = "Small cod")
#> mutate: new variable 'species2' (character) with one unique value and 0% NA

lc_margin_sad <- predict(l_cod_sad2, newdata = nd_fle, se_fit = TRUE, re_form = NA, re_form_iid = NA) %>% 
  mutate(species2 = "Large cod")
#> mutate: new variable 'species2' (character) with one unique value and 0% NA

f_margin_sad <- predict(fle_sad2, newdata = nd_fle, se_fit = TRUE, re_form = NA, re_form_iid = NA) %>% 
  mutate(species2 = "Flounder")
#> mutate: new variable 'species2' (character) with one unique value and 0% NA

cond_effects <- bind_rows(sc_margin_sad, lc_margin_sad, f_margin_sad)

ggplot(cond_effects, aes(density_flounder, exp(est),
                         ymin = exp(est - 1.96 * est_se),
                         ymax = exp(est + 1.96 * est_se), 
                         color = species2,
                         fill = species2,
                         linetype = species2)) +
  geom_ribbon(alpha = 0.2, color = NA) +
  theme(aspect.ratio = 1) +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_linetype_manual(values = c(3, 2, 1)) +
  guides(linetype = "none",
         color = guide_legend(override.aes = list(linetype = c(3, 2, 1)))) + 
  geom_line(size = 1.3) +
  labs(x = "Flounder biomass density",
       y = "Saduria feeding ratio") + 
  theme(legend.position = c(0.9, 0.9))

  
ggsave("figures/flounder_saduria_conditional.pdf", width = 17, height = 17, units = "cm")